docs/rxd-tutorials/combining currents from mod files with rxd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: combining MOD file ion channels with rxd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A version of this notebook may be run online via Google Colab at https://tinyurl.com/neuron-rxd-and-mod\n",
    " (make a copy or open in playground mode)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NEURON's reaction-diffusion infrastructure can be used to readily allow intracellular concentrations to respond to currents generated in MOD files, as long as:\n",
    "\n",
    "- `nrn_region='i'` is specified for the `rxd.Region` (so that it knows it corresponds to the electrophysiology region of the inside of the cell), AND\n",
    "the name and charge of the ion/etc are given in the `rxd.Species` declaration.\n",
    "Satisfying the above two rules also allows MOD files to see intracellular concentrations.\n",
    "\n",
    "- 3D extracellular concentrations also interoperate with electrophysiology automatically as long as name and charge are specified.\n",
    "\n",
    "As a simple example, we consider a model with just a single point soma, of length and diameter 10 microns, with Hodgkin-Huxley kinetics (which uses the built in mod file `hh.mod`), and dynamic sodium (declared using rxd but without any additional kinetics)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup NEURON library and imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's import our usual NEURON libraries and definitions. Remember you can use either `um` or `µm` for micron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from neuron import h, rxd\n",
    "from neuron.units import mV, ms, um, mM\n",
    "\n",
    "## needed for standard run system\n",
    "h.load_file(\"stdrun.hoc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now import `plotly`, a graphics library. (You could easily modify this code to use other graphics libraries like `matplotlib`, `plotnine`, or `bokeh`.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.graph_objects as go\n",
    "from plotly.subplots import make_subplots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## define morphology\n",
    "soma = h.Section(name=\"soma\")\n",
    "soma.L = soma.diam = 10 * um\n",
    "\n",
    "## add ion channels (h.hh is built in, so always available)\n",
    "h.hh.insert(soma)\n",
    "\n",
    "## define cytosol. MUST specify nrn_region for concentrations to update\n",
    "cyt = rxd.Region([soma], name=\"cyt\", nrn_region=\"i\")\n",
    "\n",
    "## define sodium. MUST specify name and charge for concentrations to update\n",
    "na = rxd.Species(cyt, name=\"na\", charge=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we could have written `h.hh.insert(soma.wholetree())` to put Hodgkin-Huxley channels everywhere in the cell that the soma is part of, but since there is only one section, this is not required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also add an excitatory synapse to receive events (these will trigger the cell to spike)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "syn = h.ExpSyn(soma(0.5))\n",
    "syn.tau = 1 * ms\n",
    "syn.e = 0 * mV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add a stimulus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spike events themselves (two events, 15 ms apart starting at h.t=10*ms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stim = h.NetStim()\n",
    "stim.interval = 15 * ms\n",
    "stim.number = 2\n",
    "stim.start = 10 * ms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Send those events to our synapse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc = h.NetCon(stim, syn)\n",
    "nc.weight[0] = 0.001"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup recording variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = h.Vector().record(h._ref_t)\n",
    "v = h.Vector().record(soma(0.5)._ref_v)\n",
    "na_vec = h.Vector().record(soma(0.5)._ref_nai)\n",
    "ina = h.Vector().record(soma(0.5)._ref_ina)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize and run the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h.finitialize(-65 * mV)\n",
    "h.continuerun(50 * ms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = make_subplots(rows=3, cols=1)\n",
    "fig.add_trace(go.Scatter(x=t, y=v, name=\"v\"), row=1, col=1)\n",
    "fig.update_yaxes(title_text=\"v (mV)\", row=1, col=1)\n",
    "fig.add_trace(go.Scatter(x=t, y=ina, name=\"ina\"), row=2, col=1)\n",
    "fig.update_yaxes(title_text=\"ina (mA/cm^2)\", row=2, col=1)\n",
    "fig.add_trace(go.Scatter(x=t, y=na_vec, name=\"[Na+]\"), row=3, col=1)\n",
    "fig.update_xaxes(title_text=\"t (ms)\", row=3, col=1)\n",
    "fig.update_yaxes(title_text=\"[Na+] (mM)\", row=3, col=1)\n",
    "\n",
    "fig.show(renderer=\"notebook_connected\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Note: On homeostatic mechanisms or the lack thereof\n",
    "\n",
    "Without any additional homeostatic mechanisms (Hodgkin and Huxley did not model sodium concentration so they did not need to include homeostatic mechanisms for it), intracellular sodium concentration will not return to baseline, and each spike will move intracellular sodium concentration closer to the extracellular concentration. Potassium concentration in this model is constant as we did not declare a potassium rxd.Species, but if we did it would also approach its extracellular concentration with each spike, until eventually the cell is not able to fire action potentials anymore."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finally: The mod file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For any distributed mechanism, including the built-in `h.hh`, one can always get the source code for the corresponding mod file using the `.code` property; e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(h.hh.code)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}